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A theoretical study of the magnetoelectronic properties of zigzag and armchair bilayer graphene 
nanoribbons (BGNs) is presented. Using the recursive Green's function method, we study the 
band structure of BGNs in uniform perpendicular magnetic fields and discuss the zero-temperature 
conductance for the corresponding clean systems. The conductance quantized as 2{n + l)Go for the 
zigzag edges and nGo for the armchair edges with Go = 2e^ /h being the conductance unit and n an 
integer. Special attention is paid to the effects of edge disorder. As in the case of monolayer graphene 
nanoribbons (GNR), a small degree of edge disorder is already sufficient to induce a transport gap 
around the neutrality point. We further perform comparative studies of the transport gap Eg and 
the localization length ^ in bilayer and monolayer nanoribbons. While for the GNRs E^'^^ ~ ^/W, 
the corresponding transport gap E^^'^ for the bilayer ribbons shows a more rapid decrease as the 
ribbon width W is increased. We also demonstrate that the evolution of localization lengths with 
the Fermi energy shows two distinct regimes. Inside the transport gap, ^ is essentially independent 
on energy and the states in the BGNs are significantly less localized than those in the corresponding 
GNRs. Outside the transport gap ^ grows rapidly as the Fermi energy increases and becomes very 
similar for BGNs and GNRs. 



PACS numbers: 81.05.Uw, 73.23.-b, 73.21.Ac 

I. INTRODUCTION 

The recent successful fabrication of monolayer 
graphenei has ignited tremendous interest because it not 
only represents a platform to model relativistic particles 
in a condensed matter material, but also provides a po- 
tential building block for future nanoelectronics. The 
honeycomb lattice of the graphene sheet imposes a lin- 
ear low-energy electronic spectrum and the correspond- 
ing extraordinary electronic behavior of the excitations 
near the Dirac point, namely massless Dirac fermions^ii^. 
Various aspects of graphene like band structures^ and 
resulting transport properties^i^i^iiiS have been studied. 
Electrostatic potentials^ and magnetic barrier a^^'^^i^^ 
have been suggested as ways to achieve tunable confine- 
ment. Moreover, there have been suggestions to induce a 
bandgap and manipulate efficiently the resistance by, for 
example, chemical dopingi^, edge disorder—, while stud- 
ies of the magnetoelectronic properties also revealed the 
anomalous integer and fractional quantum Hall effects in 
graphenei^^iia. 

More recently, increasing attention has been paid 
to graphene multilayers, in particular to bilayer 
systems2£ii^i22ii^i2i which consist of two coupled 
graphene sheets with two sublattices, A and B, in each 
layer. They are typically stacked in the Bernal form, 
where A' sites belonging to the upper layer are located 
exactly on the top of B sites to the lower layer, and B' or 
A sites are above or below the center of hexagons in the 
other layer, as shown in Fig. [TJ Bilayers show quite dif- 
ferent properties than monolayers in many respects, such 



as the quantum Hall effec t^^'^° , edge states^^, and weak 
localization^^. Graphene bilayers are also inevitably in- 
fluenced by the disorder present in the environment. The 
role of disorder in graphene bilayers has been studied 
theoretically^! and the minimal conductivity has been 
addressed by different author a^^i^^ . An important prop- 
erty of graphene bilayers is that an application of an 
electric field between the layers allows to open up and 
tune a band gap, which is highly relevant for the pos- 
sible applications in nanoelectronics. Experimentally, 
a double-gate configuration can impose a perpendicular 
electric field onto a bilayer, and a controlled transition 
from a zero-gap semiconductor to an insulator has been 
observed, which provides a direct evidence of opening a 
bias-induced band gap in a bilayer—. ft has been shown 
that a bias voltage can continuously tune a bandgap 
from zero to mid-infrared frequencies around the zero- 
energy poin1^, modify the charge density distribution in 
the graphene planes^^ and even induce the low-density 
ferromagnetism^. 

To study the properties of bulk bilayer graphene theo- 
retically, the effective Dirac-like Hamiltonian as the con- 
tinuum limit of the tight-binding model close to the Dirac 
points was mostly used, ft has been applied to the 
study of the energy spectrum of bilayers^, to the trans- 
port phenomenology through biased^ and unbiased^ 
graphene bilayers with bulk disorder as well as to the 
bilayer/monolayer interface*^. However, the energy spec- 
tra and the transport properties of bilayer graphene 
nanoribbons (BGNs), which form the focus of this pa- 
per, have not yet been addressed in the literature to the 
best of our knowledge. 
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The purpose of the present work is twofold. First, we 
develop a computational method and present a theoret- 
ical study of the magneto-band structure, Bloch states 
and conductance quantization in zigzag and armchair 
BGNs. Our method is based on the recursive Green's 
function technique that we recently developed to calcu- 
late the electronic structure and conductance of graphcnc 
nanoribbonsi^. The special feature of this technique is 
that, in contrast to earlier Green's function methods, it 
does not require self-consistent calculations for the sur- 
face Green's function, which is instead computed as a 
solution of an eigenequation. This greatly reduces the 
computation time. Making use of this technique, one 
can study the band structure and conductance of BGNs 
of various geometries and edge terminations in various 
transport scenarios including a perpendicular magnetic 
field, a bias voltage, etc. It should also be noted that 
the electron-electron interaction and screening can signif- 
icantly affect and modify electronic properties of single- 
layer graphene nanoribbon a'^^i'^^ . It is expected that 
screening in the BGNs differs from that one in GNRs, 
because the electrons in each layer are affected not only 
by the electrons belonging to the same layer, but also by 
the electrons in the second layer. A detailed knowledge 
of the electronic structure of an ideal BGNs (without 
interaction) is a prerequisite for any studies of interac- 
tion effects in this system. While accounting for electron 
interactions will be the subject of future studies, knowl- 
edge of electronic structure of the ideal BGN represents 
a necessary first step in this direction. 

Second, the effects of edge disorder on the conductance 
of BGNs are investigated. It is well established that the 
edge disorder strongly affects the transport properties of 
monolayer graphene nanoribbons leading to the suppres- 
sion of the conductance quantization and opening of the 
transport gapS^OSdOiii^^^^ie^ ^^^^^ of 

any work addressing the effect of edge disorder on the 
transport properties of BGNs and our study represents a 
first step in this direction. We demonstrate that conduc- 
tance of realistic BGNs with edge disorder share many 
similar features with that one of the monolayer GNRs. 
There are however important and interesting differences 
that we discuss in detail. 

The paper is organized as follows. In Section II, we 
sketch the geometry of the devices and the Green's func- 
tions technique that we use. In Section III, we present 
the band structures and the conductance of ideal BGNs. 
Section IV focuses on the effects of edge disorder. A 
summary and conclusions constitute Section V. 





FIG. 1; (Color online) Schematic geometry of the two- 
terminal structure for the case of (a) zigzag and (b) armchair 
BGNs. The devices under consideration are indicated by the 
shaded regions which connect to the two semi-infinite leads 
opening to the left and right. Bonds in the bottom layer 
are indicated by dashed lines and in the top layer by solid 
lines. Empty and filled circles denote the atoms in the lower 
and upper layers, respectively. Unit cells of the BGNs are 
marked by bold dashed rectangles. Vertical dashed lines with 
numbers label the slices of a unit cell; the numbers with and 
without prime labels refer to respective upper and lower lay- 
ers, for example, 1 and 1' correspond to the first slice, a is 
the carbon-carbon bond length. Numbers in the vertical di- 
rection label the sites in a slice, (c) 3D view of the graphene 
bilayer with all coupling energies. Two interpenetrating trian- 
gular sublattices sublattices in the lower and upper layers are 
denoted A, B, and A' , B' , respectively. 70 is the coupling be- 
tween nearest-neighboring sites in the same layer. Interlayer 
coupling energies include 71 between A' and B, 73 between 
B' and A, and 74 between A'{B') and A{B) (74 is neglected 
in our calculations). 



II. MODEL DESCRIPTIONS 

We shall consider the two-terminal BGN structure 
sketched in Fig. [ija) and (b). The semi-infinite left 
and right leads are made of zigzag or armchair BGNs. In 
the shaded central region, one can define a structure of 
interest, such as an electrostatic barrier, point defects, or 



edge roughness. In our approach we describe the bilayer 
by the tight-binding Hamiltonian^ 

H = ^ [Vi^iaj^^ai^i + Vtjbjjbgj) - -/Q ^ (aj.bij+h.c.) 
- li^{altb2,r + h.c.)--/3j2iKi^2,j+h.c.) (1) 
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where a/^ (b^i) is the creation operator at sublattice A 
(B), in the layer £ = 1,2, at site Ri. Vi^i is the on- 
site electrostatic potential. We use the common graphite 
nomenclature for the coupling parameters^ (see illustra- 
tion in Fig. [T] (c)): 70 = 3.16 eV is the intralayer 
nearest-neighbor coupling energy, 71 = 0.39 cV is the 
coupling energy between sublattice B and A' in different 
graphcnc layer, and 73 = 0.315 eV the hopping energy 
between sublattice A and B' in the lower and upper lay- 
ers, respectively. The other coupling energy between the 
nearest-neighboring layers, 74 w 0.04 eV, is very small 
compared with 70 and ignored below. In the presence 
of an external perpendicular magnetic field B, the hop- 
ping integral acquires the Peierls phase factor such that 
the coupling 70(3) is modified to 70(3) exp (ieS^.r' /^) with 

dr,r' = /J^ A • dl and 71 leaves unchanged, where the 
Landau gauge, A = {—By, 0, 0) was used. 

In the Green's function method, it is convenient to 
write the Hamiltonian H in the form 

H = Y,[h,] + U, (2) 

i 

where hi describes the Hamiltonian of the i-th slice, and 
U describes hopping between all neighboring slices. Each 
i-th slice of the BGN consists of two slices of the length 
N belonging respectively to the upper and the lower lay- 
ers, such that the dimension of the matrixes hi and U 
is 2N X 2N. The choices of slices are indicated in the 
Figs. [BJa) and[IJb). The numbers label the indices of 
slices in the lower layer, while the numbers with primes 
label the indices of slices in the upper layer. In the same 
slice of zigzag BGNs, the transverse sites of the upper 
layer are sitting above the sites of the lower layer. For 
armchair BGNs, we choose the slices in which the trans- 
verse sites of the lower layer are shifted by a distance of 
a/2 with respect to the sites belonging to the upper layer 
in order to include all the coupling between the nearest- 
neighboring slices. The Hamiltonian matrix of the i-tli 
slice of the BGN consists of two N x N sub-matrices de- 
scribing the slices of the upper and lower layers, which 
constitute the diagonal blocks. The off-diagonal blocks 
of the Hamiltonian matrix account for the interactions 
between the slices of the upper and lower layers. The di- 
agonal blocks of the hopping matrixes U describe the in- 
tralayer hopping between the nearest-neighboring slices, 
while the off-diagonal blocks describe the hopping be- 
tween the nearest-neighboring slices in the different lay- 
ers. Explicit expressions for the matrices [hi] and U are 
given in the Appendix TK] 

To calculate the band structures, we first study an infi- 
nite ideal BGN by using the Green's function formalism. 
A unit cell of the BGN structure consists of M slices, 
with M ^ 2 for the zigzag and A/ = 4 for the arm- 
chair as shown in Figs. [T] (a),(b). The Bloch states and 
the wave functions of BGNs are obtained by solving an 
eigenvalue problem which is formulated in the Appendix 
IB] The wave functions of the slices arc further used to 



calculate the surface Green's functions accounting for the 
effects of semi-infinite leads. 

In order to calculate the transmission coefficient we di- 
vide the structure into three regions. The left and right 
leads are modeled by ideal semi-infinite BGNs, which dif- 
fer only in that one extends to the left and the other to 
the right. These two leads are then connected to one 
other via the central device which is indicated by the 
shaded region in Figs. [1] (a),(b) (where scattering oc- 
curs). The transmission and reflection amplitudes are 
related to the total Green's functions which are com- 
puted using the standard recursive technique based on 
the Dyson equatio n^^'^" . The relevant equations are for- 
mulated in the Appendix |B] and their detailed deriva- 
tions (for the case of a monolayer graphene nanoribbons) 
can be found in Ref. p^ . The zero-temperature con- 
ductance is then calculated using the Landaucr-Biittiker 
formalism, 

where {T)pa is the transmission from incoming state a 
in the left lead to outgoing state /3 in the right lead. 

III. MAGNETO-BAND STRUCTURE AND 
CONDUCTANCE QUANTIZATION OF IDEAL 
BILAYER NANORIBBONS 

A. Magneto-band structure 

In this section we calculate the band structures of 
BGNs for different edges as well as in perpendicular mag- 
netic fields, characterized by the magnetic fiux (/) through 
a unit cell of the graphene lattice in units of the fiux 
quantum (po = h/e. 

In Figs. and [3] we present the energy band struc- 
tures of respectively zigzag and armchair BGNs with N ^ 
20, 19 and the corresponding wave functions for the val- 
ues of magnetic flux 0/0o = 0, 1/200, 1/100, 1/50. Dis- 
persion relations of the corresponding monolayer GNRs 
are plotted for comparison. As expected, the band 
structure of BGNs corresponds to two superimposed and 
somewhat deformed band structures of individual mono- 
layer GNRs. To understand this character, we first as- 
sume that the two graphene layers of a bilayer system are 
decoupled. The system is therefore expected to behave 
as two isolated monolayer GNRs with all bands in the 
dispersion relation being doubly degenerate. By turning 
on the coupling between the two layers the degeneracy of 
the states is then lifted by the presence of the coupling, 
with the degree of splitting depending on the coupling 
strength. This is in a close analogy to a textbook exam- 
ple of a coupled quantum wells, where the splitting of the 
quantum well states depends on the strength of the bar- 
riers separating the wells. For the case of BGNs the cou- 
pling between layers is given by the hopping integrals 71 



Monolayer (N=20) Bilayer (N=20) 






k(l/a) 



10 15 
Sites 



(a) 



3 
2 

T 

^ 

w -1 

-2 
-3 

(b) 3 

2 

^ i 



-1 

-2 
-3 

3 
2 
_ 1 

'o 

-1 

_2 

-3 

3 
2 

1 

-0 
-1 
_2 

-3 



Monolayer (N=19) Bilayer CN=19) 





vwm 



(c) 



fl- 




ed) 












Upper 



Lower 



1 -0.5 0.5 1 -1 -0.5 0.5 1 2 6 10 14 1^ 
k(l/a) k(l/a) Sites 



FIG. 2; The energy band structures of zigzag monolayer (left 
panel) and bilayer (right panel) ribbons with A'' = 20 for 
different magnetic flux through a hexagon (a) (?!)/0o = 0, (b) 
ct>/<po = 1/200, (c) 0/<j!>o = 1/100 and (d) cj>/cj>o = 1/50. The 
corresponding wave functions at the energy represented by 
full circles are shown to the right (the top flgure for each <l>/4'o 
shows the wave function of the monolayer, while below, the 
wave functions of the upper and the lower layer are shown). 



and 73 whose magnitudes (~ 0.3 — 0.4 eV) determine the 
value of splitting of the corresponding monolayer bands. 

The features of the band structure of BGNs at i? = 
outlined above persist in the presence of a magnetic field. 
As B increases to 0/00 = 1/200 the bands close to the 
Dirac points begin to flatten. This represents the onset 
of Landau level formation when the classical cyclotron 
radius rc for the low-velocity states (i.e. for those close 
to the Dirac points) becomes smaller than the ribbon 
width ( rc = hk/eB = kl%, h = y^h/eB (=1.3 nm 
at 0/00 = 1/200) being the magnetic length). As the 
field is increased to 0/0o = 1/50, the cyclotron motion 
dominates and the Landau levels become well developed. 

It is interesting to note that similar to the case of 
monolayer GNRs, partly flat bands at E = exist in 
the band structures of the zigzag BNGs corresponding 
to the localized states on each edge, see Fig. O Com- 
pared to the monolayer cases, the number of edge states 
is doubled for each edge in bilayer ribbons. A detailed 
theoretical treatment of edge states in bilayer systems 
has been made recently^. 

Figures I2I3I also show some representative wave func- 
tions for different magnetic fields. As expected, the wave 
functions in the BGNs closely follow the corresponding 



FIG. 3: The energy band structures of metallic armchair 
monolayer (left panel) and bilayer (right panel) ribbons with 
iV = 19 for (a) 0/0o = 0, (b) 0/0o = 1/200, (c) 0/0o = 1/100 
and (d) 0/00 = 1/50. The corresponding wave functions are 
shown to the right. 



wave functions of the monolayer GNRs. As the magnetic 
field increases, the wave functions are pushed towards the 
ribbon edges, signaling the formation of the familiar mag- 
netic edge states which correspond to classical skipping 
orbits^. 

It should be noted that our nanoribbons are rather 
narrow and therefore the onset of the formation of the 
Landau levels (when the magnetic length becomes 
comparable to the ribbon's width) can be reached only 
for an unrealistically high magnetic fields (for example, 
0/00 = 1/50 corresponds to i? = 400 T). In realistic 
experimental samples of submicron width, the formation 
of Landau levels scales down to much smaller magnetic 
fields. Due to computational limitations, however, it is 
rather impractical to consider wider ribbons and in the 
present study we therefore consider the same physics at 
higher fields. 

The electronic properties of armchair graphcne 
nanoribbons depend sensitively on the ribbon width. The 
GNRs arc metallic when iV — 1 is a multiple of 3; other- 
wise, there is an energy gap between the conduction and 
valence band. Fig. HJa-b) displays the band structures of 
the armchair BGN and GNR at B = with N = 8,7. It 
is evident that both the monolayer and bilayer graphene 
ribbons with N = 8 possess a bandgap around the zero- 
energy point. However, the energy gap of the BGN is 
narrower than that of the monolayer GNR even though 
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FIG. 4: The energy bands of semiconducting and metallic 
armchair nanoribbons in B = 0: (a) monolayer with N = 8; 
(b) bilayer with N = 8. In (c) and (d), the band structures 
of the monolayer and bilayer systems with N = 7 are shown. 



they have the identical width. For the armchair BGN 
and GNR with N = 7, the band structures show metal- 
lic characteristics as shown in Fig. ID^c) and (d). In 
contrast to the linear dispersion relation near the Dirac 
point in the monolayer GNR (Fig. [H^c)), the armchair 
BGN (Fig. HKd)) shows a parabolic dispersion relation 
near the zero-energy point, which leads to new elemen- 
tary excitations called massive Dirac fermions. These 
massive non-relativistic particles are chiral and can be 
described by an ofF-diagonal, Dirac-like Hamiltonian as 
the continuum limit of the tight-binding model^S. This 
Hamiltonian of the bulk bilayer system gives rise to two 
parallel parabolic bands separated by an amount of ~ 71 
close to the zero-energy point, which is consistent with 
our observations for nanoribbon structures. 



B. Conductance quantization 
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FIG. 5: Upper panels: the energy spectra for the zigzag (a) 
and armchair (c) bilayer ribbons with A'^ = 46. (For the zigzag 
ribbon only the branches in the vicinity of one Dirac point are 
shown). Lower panels: The energy dependence of the zero- 
temperature conductance of the ideal zigzag (b) and armchair 
(d) ribbons associated with the energy spectra. 



Our calculations shows that the corresponding conduc- 
tance of the zigzag and armchair bilayer ribbons is given 
respectively by Cf — 2{n+ l)Go and G^qj^ = nGo, 
see Fig. [5] The minimum conductance of an ideal zigzag 
BGN is 2Go, whereas the minimum conductance of an 
ideal metallic armchair BGN is Gq. 

Note that the conductance steps are not equidistant 
along the energy axis. For example, for the armchair 
BGN steps at Go, 3Go, and 5Go are more pronounced 
than those at 2Go, 4Go, and 6G0. The widths of conduc- 
tance steps are determined by the energy scale between 
the successive modes in the energy spectrum, which in 
turn depends on the ribbon width and the energy inter- 
val. In wider armchair BGNs, some conductance quanti- 
zation steps are poorly resolved or even unresolved which 
might lead to an apparent quantization in units of 2Go. 



IV. EFFECT OF EDGE DISORDER ON THE 
CONDUCTANCE AND LOCALIZATION 
LENGTH 



Now we are in the position to study the transport prop- 
erties of BGNs. In clean ideal nanoribbons the conduc- 
tance is determined by the number of propagating modes 
at the Fermi energy. Each spin-degenerate propagating 
mode contributes to the total conductance by the con- 
ductance unit Go = 2e? /h. Figure O shows the energy 
spectra for zigzag and armchair BGNs with = 46 and 
a corresponding conductance as a function of the Fermi 
energy at _B = 0. Note that analytical expressions for the 
conductance quantization in monolayer GNR were pro- 
vided by Onipko^ who showed that the conductance of 
the monolayer GNRs is given by G§^^ = 2(n + l/2)Go 
for the zigzag ribbons and Gqj^j^ = nGo for the arm- 
chair ribbons, where the integer number n = 0,1,2,---. 



In this section we calculate the conductance of real- 
istic BGNs with edge defects that are created by ran- 
dom removal of carbon atoms from the edges of the up- 
per and lower layers. It was shown previously that the 
edge disorder strongly affects the conductance of mono- 
layer GNRs^^^^^^^^^i^^i^ . In particular, it 
has been demonstrated that even very modest edge dis- 
order is sufficient to strongly suppress the conductance 
and to induce the conduction energy gap in the vicin- 
ity of the Dirac point and to lift any difference in the 
conductance between nanoribbons of different edge ge- 
ometry (i.e. zigzag and armchair) ^^1^^ . This was related 
to the pronounced edge-disorder-induced Anderson-type 
localization which leads to a strongly enhanced density 



6 



L= 25 nm 



Zigzag 



L= 123 nm 



■5 O 



14 
10 
6 

2 

30 



AO ^^^f 

— Ai r 

A2 • ■' 

A3 != . 

J ^ 









L = 22 nm Armchair 




(a) bilayer 



(b) monolayer 



L= 107 nm 





^1 



































0.1 0.2 0.3 0.4 0.5 0.1 0.2 0.3 0.4 0.5 

Ef(To) Ep(yo) 

FIG. 6: (Color online) The conductance as a function of the 
Fermi energy with varying edge disorder for different widths 
and lengths. Three typical disorder strengths are considered: 
AO: perfect edges; AI: pi = 0.1, p2 = Ps — 0; A2: pi = 0.3, 
P2 = 0.2 and ps = 0.1; A3: pi = 0.5, P2 = 0.4 and pa = 0.3. 
The insets show the positions of the atoms for removing with 
the probabilities pi,P2, and pa. (a)-(d) the zigzag BGNs with 
widths 5 nm-10 nm and length 25 nm-123 nm. (e)-(h) the 
armchair BGNs with widths 5 nm-11 nm and lengths 22 nm- 
107 nm. 



of states at the edges and to blocking of conductive paths 
through the ribbonsi^. In the present study we use the 
model for the edge disorder applied previously to mono- 
layer GNRsi^ and compare the effect of the edge disorder 
on the conductance and the localization length of mono- 
layer GNRs and BGNs. We model the missing atoms 
by setting the corresponding hopping matrix elements to 
zero. The edge roughness is controlled by three probabil- 
ities pi, p2, P3- Pi is the probability of a missing atom in 
the outermost row; pi(i = 2, 3) is the conditional proba- 
bility for a missing atom in the z-th row away from the 
edge if at least one adjacent atom in row — 1) is missing, 
see the illustration in the inset of Fig. [HI As graphene 
is known to have few bulk defects in general^ we do not 
consider bulk vacancies. We also disregard the effect of 
hydrogen capture by dangling bonds at the edge, since 
it has been shown to be of minor importance for ribbons 
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FIG. 7: The inverse of the energy gap versus ribbon width for 
(a) armchair BGNs and (b) armchair monolayer GNRs. The 
energy gap Eg is defined as the interval where G < 10""^ x 
2eVft. 



wider than a few nanometers^i^. 

Figure [6] shows the conductance of zigzag and arm- 
chair BGNs of different lengths and widths as a func- 
tion of the Fermi energy for four edge disorder strengths. 
The conductance of the BGNs shows the same qualitative 
features as the conductance of the monolayer GNR&i^. 
The main feature is that a relatively small edge disor- 
der strongly suppresses the conductance and completely 
destroys the quantization steps for both zigzag and arm- 
chair BGNs. As a result, no difference in the conductance 
is expected for realistic zigzag and armchair BGNs. Arm- 
chair BGNs show a well-pronounced transport gap in the 
vicinity of the charge neutrality point. As in the case 
of monolayer ribbons, the zigzag BGNs are more robust 
to this disorder-induced suppression of the conductance. 
This behavior has been related to the presence of the 
edge states in the zigzag nanoribbon close to the Dirac 
point^. 

Recent experiments and calculations have shown that 
the transport gap in monolayer GNRs is approxi- 
mately proportional to the inverse of the ribbon width 
^_9 ^i5,4i^44,46 ^ In Fig. [7] wc plot the inverse of the trans- 
port gap, E~^, for the armchair BGNs and monolayer 
GNRs as a function of the ribbon width W. While for 
the monolayer ribbons E^^^ ~ ^/W, the correspond- 
ing transport gap E^'^^ for the bilayer ribbons shows a 
more rapid decrease as W increases. This indicates that 
edge-disorder induced locaHzation of the states is less pro- 
nounced in BGNs in comparison to GNRs. In order to 
elucidate this point further, we have performed a com- 
parative study of the localization length ^ as a function of 
the Fermi energy Epm both systems, choosing nanorib- 
bons of width W ^ 11 nm, corresponding to 46 atomic 
sites in y-direction. The vacancy probabilities have been 
set to pi = 0.5, p2 = 0.4 and p^ — 0.3, respectively. 
The nanoribbons are thus identical to that one of Fig. [6] 
(h), configuration A3. S,{Ef) is determined by calculat- 
ing the conductance G as a function of the nanoribbon 
length L for an ensemble of 5000 disorder configurations 
for each length. As the energy increases, the system is 
expected to undergo a transition from the strongly local- 
ized regime where the conductance decays exponentially 
as L increases, to the Ohmic regime where G oc 1/ L. The 
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conductance is expected to obey the scaling law 
ln(l + 1/g) = L/^, 



(4) 



with the dimensionless conductance g ~ Gjie^jnK^. 
We find that for all energies, the simulated conductances 
can be fitted with high accuracy by Eq. as exem- 
plified in the inset of Fig. EIc), with ^ being the fit 
parameter. Fig. [5] (c) shows the localization length de- 
termined that way as a function of the Fermi energy for 
both monolayer and bilayer ribbons. Two regimes can 
be identified, depending on whether the Fermi energy 
Ep is inside or outside the transport gap. For energies 
inside the transport gap, \Ep\ < Eg, (with Eg « O.O670 
for the particular system under study), ^ is essentially 
independent of Ep and comparable to W. In this inter- 
val, ^ in the BGN is roughly a factor of 2 larger than in 
the GNR. As Ep increases above the transport gap, ^ in- 
creases rapidly, while the difference in ^ for the mono-and 
the bilayer system fluctuates around zero. We emphasize 
that these fluctuations are well above the noise level, due 
to the large ensemble used. 

This behavior can be interpreted qualitatively as fol- 
lows. Inside the transport gap, \Ep\ < Eg, large disorder- 
induced potential barriers are present in each layer which 
hamper the electron transfer along the wire. If an elec- 
tron encounters such a barrier in one layer, it can bypass 
it by hopping to the second layer. Because such a bypass- 
ing is apparently not possible in a system with a single 
layer the electrons are less localized in the BGNs. 

Outside the transport gap, \Ep\ > Eg, conducting 
channels open up in both layers and therefore the in- 
terlayer transfer becomes less likely compared to the in- 
tralayer transfers. Because of this, electronic transport 
in each layer is much less affected by the presence of 
the neighboring level which results in similar localization 
lengths for the monolayer and the bilayer nanoribbons. 
The differences in ^ for the mono- and the bilayer sys- 
tem mentioned above can be traced to the features in 
the dispersion relations corresponding to the opening of 
new propagating channels. A correlation between the 
dispersion relation and the behavior of S,{E) is clearly 
seen in Fig. [SI For example, the strong increase of ^ be- 
tween O.O670 < E < O.I70 of the monolayer is related to 
the occupation of the second mode. As the energy (and 
therefore a number of modes) increases, the difference of 
$,{E) between mono- and bilayer nanoribbons diminishes. 

An inspection of a typical local density of states 
(LDOS) pattern of a representative member of the 
nanoribbon ensemble substantiates this interpretation of 
the properties of £,{E), see Fig. [HI We first generate two 
monolayer GNRs with the same concentration but differ- 
ent configurations of edge disorder, and then couple them 
to form the lower and the upper layers of a bilayer BGN. 
We then separately calculate the LDOS for these three 
structures (i.e. for two different monolayer GNRs and 
one BGR). The pronounced enhancement of the LDOS 
at the edges (note the logarithmic scale) and the sub- 
stantial fluctuations at larger scales indicate Anderson 
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FIG. 8: Bottom: The localization length ^ as a function of 
the Fermi energy for a ribbon of the width W = 11 nm. Inset: 
length dependence of ln(l -I- 1/g). Top: The monolayer and 
bilayer dispersion relations for the corresponding nanoribbons 
with perfect edges. 



Ep=0.05T„ 
(a) Monolayer #1 









(b) Lower layer 














(c) Upper layer 


r ' 






p 






(d) Monolayer #2 





20 

X (nm) 



E„-0.09T„ 
(e)Monolayer#l 



(f) Lower layer 



(g) Upper layer 



(h) Monolayer #2 



10 20 30 40 

X (nm) 



FIG. 9: (Color online) The local density of states in repre- 
sentative parts of the edge disorder regions, (a) and (d) are 
the LDOS of two independent monolayer at Ep = O.O570 
(i.e. inside the transport energy gap) and evolve into (b) and 
(c), respectively, when the interlayer interaction is turned on. 
(e) and (h) are the LDOS of two independent monolayer at 
Ef = 0.0970 (i.e. outside the energy gap) and evolve into (f) 
and (g), respectively, when the interlayer interaction is turned 
on. 
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localization as discussed in detail clsewherei^. 

Figures [9] (a)- (d) show the LDOS plots for Ep lying in- 
side the energy gap. A comparison of the LDOS pattern 
shows an anticorrelation between the lower and upper 
layers of the BGN (i.e. the enhanced LDOS in one layer 
is accompanied by the suppressed LDOS in the second 
layer). At the same time, there is not much correlation 
between LDOS patterns of the BGN and the correspond- 
ing LDOSs of the monolayer GNRs. This is consistent 
with the interpretation presented above where the inter- 
layer transfer is dominant such that the electrons en- 
countering a potential barrier (which reduces the LDOS) 
in one layer, jump to the second layer and enhance the 
LDOS of the neighboring sites therein. 

The LDOS patterns for a representative value of Ep 
above the transport gap demonstrate that the intralayer 
transfer is dominant in this regime. Indeed, in this case 
a correlation can be seen between the LDOS patterns of 
the BGN and the corresponding LDOSs of the monolayer 
GNRs which signifies the presence of the open propagat- 
ing channels in each layer. Note that there is also a cer- 
tain correlation even between two layers in the BGN. For 
example, the lower layer of the BGN (Fig. Eljf)) shows 
some additional LDOS maxima which are not present in 
its constitutive monolayer #1 (see Fig. [HJe)) but cor- 
relate to the LDOS maxima in the second layer of the 
BGN (Fig. [9] (g)). This correlation between layers in- 
dicates that due to the coupling between the layers the 
defects in the upper layer simply generate additional de- 
fects of comparable strength in the lower layer and vice 
versa. 

Finally we note that no attempts were made to quan- 
tify the correlations in the LDOS patterns. We rather 
restricted ourselves to a visual inspection which explain 
qualitatively the behavior of £,{Ef). 



V. SUMMARY AND CONCLUSIONS 



disorder strongly suppress the conductance of the BGN, 
introduces the transport energy gap in the vicinity of 
the charge neutrality point and completely destroys the 
quantization steps for both zigzag and armchair BGNs. 
We calculate the localization length as a function of the 
Fermi energy Ep and identify two distinct regimes de- 
pending on whether Ep is inside or outside the transport 
gap. The localization length inside the transport gap is 
larger in BGNs than in GNRs with identical edge rough- 
ness. This difference, however, disappears at energies 
above the transport gap. 
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APPENDIX A: HAMILTONIAN AND COUPLING 
MATRICES 

In this appendix we provide explicit forms of the hamil- 
tonians for the i-th slices hi and coupling matrices U . 
denotes the potential on the site i in the upper layer and 
Vi the potential on i-th site in the lower layer. The diag- 
onal blocks of the matrices correspond to the Hamiltoni- 
ans of layers, and off-diagonal blocks give the intcrlayer 
interactions. 



We have studied the electronic structures and trans- 
port properties of Bernal-type graphene bilayer nanorib- 
bons using the recursive Green's function technique. The 
band structures of the zigzag and armchair BGNs with 
uniform perpendicular magnetic fields have been com- 
puted for the entire energy regime of the 7r-bands. The 
calculated band structure corresponds to two superim- 
posed and somewhat deformed band structures of indi- 
vidual monolayer GNRs. The splitting between the cor- 
responding monolayer bands is of the order of 0.3 — 0.4 
eV and is determined by the intcrlayer hopping integrals 
7i and 73 providing coupling between the layers. The 
conductance of the clean (disorder-free) zigzag system 
is quantized as 2(n -\- l)Go, whereas the clean armchair 
system is quantized as tiGq with Go = 2e^//i being the 
conductance unit and n an integer number. Further- 
more, the effect of edge disorder, which is highly rele- 
vant in realistic samples, has been investigated. As in 
the case of monolayer ribbons, a relatively small edge 
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2. Armchair edge 
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APPENDIX B: FORMALISM 



In this Appendix we provide main formulas used for 
calculation of the dispersion relations, Bloch states and 
their velocities and the transmission and reflection am- 
plitudes for the scattering problem for bilayer graphenc 
ribbons. All these formulas represent a straighforward 
generalization of the corresponding formular derived in 
Ref. [l^ for monolayer ribbons. 

The band structure is computed by solving a eigen- 
value of the form 




ikM 




(Bl) 



with M being the periodicity [M ~ 2 and 4 for the zigzag 
and armchair BGRs respectively), 



Starting from the Schrodingcr equation and using a 
definition of the group velocity, we obtain the group ve- 
locity for graphene nanoribbons 



1 X - O 



mm) 



(B2) 



M ^ dk I \ip^\ 
where the summation is performed over all slices of the 
unit cell, and 



(B3) 



is a vector composed of the matrix elements ipi^j = 
{Qaij\ip) (Note that vectors ^pi can be obtained from ipi 
via the relation ijji — e''^*(^i). 

To account for the effects of the leads, one needs to cal- 
culate the surface Green's function. The latter equations 
can be used for determination of Tr , 



1,0 



*i*o 



(B4) 



where ^'i and ^'o are the square matrixes composed of the 
matrix-columns i/'i' and , (1 < a < 2A^), Eq. (jBl[) . i.e. 
*i = (V-i,...,^?"^); *o = The expression 

for the left surface Greens function F; (i.e. the surface 
function of the semiinfinitc ribbon open to the right) is 
derived in a similar fashion. 



(B5) 



where the matrixes "I'm and 5'a/+i are defined in a sim- 
ilar way as \E'i and \l/o above. 

The transmission and reflection amplitudes is calcu- 
lated using the expressions: 



Ti 







-I 







$iT = -G^'°(J7°'i«'i - r;"^$o)- 



(B6) 



/ being the unitary matrix. 

This eigenvalue problem Eq. (jBl[) gives the eigenfunc- 
tions I 1 < a < '2N and the set of Bloch eigenvectors 
{ka} which includes both propagating and evanescent 
states. The latter can be easily identified by a non-zero 
imaginary part. In order to separate right- and left- 
propagating states, fc+ and k~ , we compute the group 
velocities of the Bloch states Va = , whose signs de- 
termine the direction of propagation ('+' stands for the 
right-propagating and '-' for the left propagating states). 



G"^"(C/"-^*i 



rr'^o) 



$0- 



(B7) 



The matrices T and R with the dimension 2A^ x A'^prop are 
composed of the transmission and reflection amplitudes 
tpa and r^Q,, (with Np^-op being the number of propagat- 
ing modes in the leads). The Green's functions G^^'^^''^ 
and G"'*^ are obtained from the standard recursive tech- 
nique based on the Dyson's equation^. The transmis- 
sion and reflection are related to their amplitudes by 

(T)i3a = vt /vi\t[3a\'^ and {R)f3a = ^'rt /^'^ k/3a P • 
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